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ANALYSIS OF COMPOSITE SKIN-STIFFENER DEBOND SPECIMENS USING 
VOLUME ELEMENTS AND A SHELL/3D MODELING TECHNIQUE 
RONALD KRUEGER AND PIERRE J. MINGUET^ 


Abstract. The debonding of a skin/stringer specimen subjected to tension was studied using three- 
dimensional volume element modeling and computational fracture mechanics. Mixed mode strain energy release 
rates were calculated from finite element results using the virtual crack closure technique. The simulations revealed 
an increase in total energy release rate in the immediate vicinity of the free edges of the specimen. Correlation of the 
computed mixed-mode strain energy release rates along the delamination front contour with a two-dimensional 
mixed-mode interlaminar fracture criterion suggested that in spite of peak total energy release rates at the free edge 
the delamination would not advance at the edges first. The qualitative prediction of the shape of the delamination 
front was confirmed by X-ray photographs of a specimen taken during testing. The good correlation between 
prediction based on analysis and experiment demonstrated the efficiency of a mixed-mode failure analysis for the 
investigation of skin/stiffener separation due to delamination in the adherents. 

The application of a shell/3D modeling technique for the simulation of skin/stringer debond in a specimen 
subjected to three-point bending is also demonstrated. The global structure was modeled with shell elements. A 
local three-dimensional model, extending to about three specimen thicknesses on either side of the delamination 
front was used to capture the details of the damaged section. Computed total strain energy release rates and mixed- 
mode ratios obtained from shell/3D simulations were in good agreement with results obtained from full solid 
models. The good correlations of the results demonstrated the effectiveness of the shell/3D modeling technique for 
the investigation of skin/stiffener separation due to delamination in the adherents. 

Key words, composite materials, delamination, finite element analysis, fracture mechanics 

Subject classification. Structures and Materials 

1. Background. Many composite components in aerospace structures are made of flat or curved panels with 
co-cured or adhesively-bonded frames and stiffeners. Recent studies which focussed on the debonding mechanism 
included testing of skin/stiffener panels and failure analysis using shell models [1-4]. A consistent step-wise 
approach has been developed over the last decade which uses experiments to detect the failure mechanism, 
computational stress analysis to determine the location of first matrix cracking and computational fracture mechanics 
to investigate the potential for delamination growth. Testing of skin gage stiffened panels designed for pressurized 
aircraft fuselage has shown that bond failure at the tip of the frame flange is an important and very likely failure 
mode [5]. Comparatively simple specimens consisting of a stringer flange bonded onto a skin have been developed 
to study skin/stiffener debonding [6-8], The failure that initiates at the tip of the flange in these specimens is 
identical to the failure observed in the full-scale panels and the frame pull-off specimens [7, 9, 10]. In a related 
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study which used the same type of simulation specimen a failure criterion based on the normal strain perpendicular 
to the fiber direction was proposed [11]. 


During previous studies two-dimensional finite element models were used because modeling and 
computational times remain affordable, especially if many different configurations have to be analyzed during the 
initial design phase [9, 10, 12]. Results from a recent study where energy release rates from two-dimensional 
analyses were compared to data from full three-dimensional simulations indicated that plane stress and plane strain 
models yield upper- and lower -bound results. For more accurate predictions, however, a three-dimensional analysis 
is required [13, 14]. For detailed modeling and analysis of the damage observed during the experiments, therefore, 
the shell/3D modeling technique appears to offer a great potential for saving modeling and computational effort 
because only a relatively small section in the vicinity of the delamination front needs to modeled with solid 
elements. The technique combines the accuracy of the full three-dimensional solution with the computational 
efficiency of a plate or shell finite element model and has been demonstrated for various applications [15-17]. In 
future applications the shell/3D modeling technique therefore is expected to be used for structural components such 
as stringer stiffened panels as outlined in the overview chart of Figure 1. 

The first objective of this study is to complement the results of an earlier investigation of the skin/stringer 
debond specimen subjected to tension, where the three-dimensional mesh in the vicinity of the free edges was not 
fine enough to accurately capture the influence of the free edges on the distribution of the energy release rates across 
the width. [13, 14]. An improved mesh with a refined zone near the free edges was generated and delaminations of 
various lengths were discretely modeled at the locations where delaminations were observed during the experiments. 
Mixed mode strain energy release rates were calculated from finite element results using the virtual crack closure 
technique [18]. Computed strain energy release rates along the delamination front contour were used to correlate the 
bonding failure in the tension test configuration using a mixed-mode interlaminar fracture criterion [19, 20]. 

The second objective of this study is to demonstrate the use of the shell/3 D modeling technique for the 
simulation of skin/stringer debond specimen subjected to three-point bending. A local three-dimensional model, 
extending to about three specimen thicknesses on either side of the delamination front was used. Delaminations of 
various lengths were discretely modeled at the locations where delaminations were observed during the experiments. 
Mixed mode strain energy release rate distributions were computed across the width of the specimens using the 
virtual crack closure technique. The results were compared to mixed mode strain energy release rates obtained from 
computations where the entire specimen had been modeled with solid elements [14]. 

2. Problem Description. The current study focuses on skin-stiffener debonding resulting from buckling of a 
thin-gage composite fuselage structure as described in reference [10]. In that study, the specimens consisted of a 
bonded skin and flange assembly as shown in Figure 2(a). An IM7/8552 graphite/epoxy system was used for both 
the skin and flange. The skin was made of prepreg tape with a measured average ply thickness of h =0. 148 mm and 
had a [45/-45/0/-45/45/90/90/-45/45/0/45/-45] lay-up. The flange was made of a plain-weave fabric with a thickness 
of h =0.212 mm. The flange lay-up was [45/0/45/0/45/0/45/0/45]f, where the subscript “f” denotes fabric, “0” 
represents a 0°-90° fabric ply and “45” represents a 0°-90° fabric ply rotated by 45°. The measured bondline 
thickness averaged 0.178mm. Specimens were 25.4-mm wide and 177.8-mm long. The properties of the 
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graphite/epoxy material and the adhesive were measured at Boeing and are part of the standard design database for 
the V-22 tilt-rotor aircraft. Typical properties are summarized in Table 1. 

Four quasi-static tension tests with a gage length of 101 .6 mm were performed as shown in Figure 2(b) [10]. 
The value of the damage onset load was averaged from four tests and determined to be P=17.8kN with a 
coefficient of variation of 8.9%. Three-point bending tests were performed with a bottom support span of 101 .6 mm 
as shown in Figure 2(c). The value of the damage onset load was averaged from four tests and determined to be 
Q =427.6 N with a coefficient of variation of 12.8%. The tests were terminated when the flange debonded from the 
skin. Damage was documented from photographs of the polished specimen edges at each of the four flange comers 
identified in Figure 2(a). Corners 1 and 4 and comers 2 and 3 showed similar damage patterns for both tests. The 
damage at corners 2 and 3, formed first and consisted of a matrix crack in the 45° skin surface ply and a 
delamination at the +45°/-45° interface as shown in Figure 3. Therefore, this damage pattern has been the focus of 
the current and related earlier analyses [10, 13]. 

Fatigue tension tests were performed at a cyclic frequency of 5 Hz, an R-ratio of 0.1 and load levels 
corresponding to 70%, 60%, 50% and 40% of the quasi-static damage onset load. The damage was monitored using 
a Questar digital microscope on one edge and an optical travelling microscope on the other edge. The specimen 
edges were painted white to make the cracks and delaminations more visible. Damage was documented based on 
location at each of the four comers identified in Figure 2(a). The number of cycles at which the first matrix crack 
appeared was recorded as well as the number of cycles to delamination onset. After the test, photographs of the 
polished specimen edges were taken under a light microscope to document the occurrence of matrix cracking and 
delamination. For two specimens, at intervals of damage growth, the test was stopped and the specimens removed 
from the grips. Photographs of the polished specimen edges were then taken under a light microscope to document 
the progression of the damage. The damage was found to be the same as observed for the quasi-static experiments 
[ 10 ]. 

The complex nature of the failure observed during the experiments, where the delamination changed across 
the specimen width from a delamination running at the skin surface 45°/-45° layer interface to a delamination 
propagating in the bondline above (see Figure 3), suggests the need for a three-dimensional model. The current 
study presents an intermediate step where three-dimensional models were created by extruding two-dimensional 
models across the width. The fact that the delamination changed across the specimen width from a delamination 
running at the skin surface 45°/-45° layer interface to a delamination propagating in the bondline above, however, is 
still not accounted for in this model. Nevertheless, the three-dimensional model takes width effects into account and 
therefore provides additional insight about the significance of edge effects. 

3, Analysis of Composite Skin-Stiffener Debond Specimens Using Volume Elements. The first 
objective of this study was to complement the results of an earlier investigation of the skin/stringer debond 
specimen subjected to tension, where the three-dimensional mesh in the vicinity of the free edges was not fine 
enough to accurately capture the influence of the free edges on the distribution of the energy release rates across the 
width. [13, 14]. 
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3.1. Three-Dimensional Finite Element Model. The deformed three-dimensional model of the specimen 
with load and boundary conditions is shown in Figure 4(a) for the tension load case. The specimen was modeled 
with ABAQUS* solid twenty-noded hexahedral elements C3D20R with quadratic shape functions and a reduced 
integration scheme. A two-dimensional mesh was made first in the x-y plane, which was then extruded across the 
width to create the model shown in Figure 4. 

A refined mesh was used in the critical area of the 45° skin ply where cracking was observed during the tests. 
Outside the refined area, the mesh was modified to prevent the three-dimensional model from becoming excessively 
large as shown in Figure 4(b). The skin plies were grouped into three layered elements with -45/45/90, 90/-45/45 
and 0/45/-45 respectively, thus taking advantage of the composite solid element option in ABAQUS® [21]. The 
fabric layers were grouped into two layered elements as shown in Figure 4(b). In the transition region at the flange 
several plies were modeled by one element with material properties smeared using the rule of mixtures [22]. This 
procedure used did not ensure the full A-B-D contribution of the plies; however, it appeared suitable for the small 
transition region to enforce a reasonable model size. 

During earlier studies a two-dimensional model had been extruded into only ten uniformly spaced elements 
across the width of the specimen [13]. The final model had 32,090 elements, 145,432 nodes yielding at total of 
436,296 degrees of freedom and required about eight hours of CPU time on a SGI Origin 3200 workstation. This 
three-dimensional mesh, however, was not fine enough in the vicinity of the free edges (z=0.0 mm and z=25.4 
mm) to accurately capture the influence of the free edges on the distribution of the energy release rates across the 
width. For the simulation of the tension test, therefore, the two-dimensional model was extruded into twenty 
elements across the width of the specimen, with a refined zone (0.74 mm, five elements) near the free edges (z=0.0 
mm and z=25.4 mm) as shown in Figures 4(b) and (c). This resulted in a model with 64,180 elements and 
286,104 nodes yielding at total of 858,3 12degrees of freedom. The nonlinear analysis required about 29 hours of 
CPU time on the same workstation. 

The Virtual Crack Closure Technique (VCCT) was used to calculate the mode I, II and III components of the 
strain energy release rates for the modeled delamination [18, 23]. For the entire investigation, the ABAQUS* 
geometric nonlinear analysis procedure was used. It was assumed that the material behaved linear elastic. This was 
done in accordance with previous studies which were used as reference solutions [14]. 

3.2. Fracture Mechanics Analysis. A fracture mechanics approach was used to investigate delamination 
onset once the initial crack had formed. During a series of nonlinear finite element analyses, strain energy release 
rates were computed at each front location for the loads applied in the experiments. A three-dimensional plot of the 
distribution of the energy release rate across the width of the specimen is shown in Figure 5. Values at the free edge 
(z=0.0 mm and z=25.4 mm) have been excluded from the plot as the model is not fine enough to accurately capture 
the energy release rates at the immediate edge. Along the length (x-coordinate) it was observed that after a small 
initial drop the computed total energy release rate increases sharply with delamination length, reaches a peak value 
and gradually decreases. Across the width (z-coordinate) the computed total energy release rate gradually increases 
with z before it first drops off near the free edges and then sharply increases in the zone of the mesh refinement in 
the immediate vicinity of the free edges (z=0.0 mm and z=25.4 mm). 
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The variation of mixed mode ratio Gs /Gy with delamination length is shown as a three-dimensional plot in 
Figure 6. Here Gs denotes the sum of the in-plane shearing components GifrGm, and G r denotes the total energy 
release rate Gri-G//+G///, where G ; is the opening mode. For two-dimensional analyses, where G n f= 0, this definition 
is equal to the previously used definition of the mixed mode ratio, Gn /Gj. For three-dimensional analysis, which 
also yields results for the scissoring mode G///, the modified definition of Gs is introduced. A three-dimensional 
plot of the distribution of the mixed-mode ratio G S /G T across the width of the specimen is shown in Figure 6. 
Values at the free edge (z=0.0 mm and z=25A mm) have been excluded from the plot as above. Along the length 
(x-coordinate) it was observed that the delamination initially starts with high shearing components, followed by a 
drop which is equivalent to an increase in opening mode I. For longer delaminations a gradual increase in shearing 
components is observed. Across the width (z-coordinate) the computed mixed mode ratio indicated high shearing 
components near the edges, followed by a drop toward the center of the specimen, which is equivalent to an increase 
in opening mode I. 

The distribution of the computed total energy release rate Gy across the width (z-coordinate) of the specimen 
is plotted in Figure 7 for the longest delamination modeled (jc= 3 1.2 mm). Across the width the computed total 
energy release rate gradually increases with z before it first drops off near the free edges and then sharply increases in 
the zone of the mesh refinement in the immediate vicinity of the free edges (z=0.0 mm and z=25.4 mm). Results 
from a previous analysis using a mesh with 10 elements uniformly spaced across the width have been included in 
the plot for comparison [14]. As expected the results obtained from the model with the refined mesh near the edges 
yield additional information about the steep increase in Gy which could not be captured with the model used earlier. 
The distribution of the mixed-mode ratio Gs /Gy across the width (z-coordinate) of the specimen is plotted in 
Figure 8 for the same location (*=31.2 mm). Across the width the computed mixed mode ratio shows high shearing 
components near the edges, followed by a drop toward the center of the specimen, which is equivalent to an increase 
in opening mode I. As before the results form an earlier study were included in the plot of Figure 8. The current 
model provided a smoother distribution across the width and yielded additional information about the steep increase 
in mixed-mode ratio Gy which could not be captured with the model used earlier. 

3.3. Mixed-Mode Failure Investigation. Delamination onset or growth for two-dimensional problems, is 
predicted by comparing calculated mixed-mode energy release rate components to interlaminar fracture toughness 
properties measured over a range of mode mixities from pure mode I loading to pure mode II loading [24-27]. A 
quasi static mixed-mode fracture criterion is determined by plotting the interlaminar fracture toughness, G r , versus 
the mixed-mode ratio, G/j/Gt, determined from data generated using pure Mode I DCB (G///Gr=0), pure Mode II 
4ENF (G//Gt= 1), and mixed-mode MMB tests of varying ratios, as shown in Figure 9 for IM7/8552 [20]. A curv e 
fit of these data is performed to determine a mathematical relationship between G r and G///Gy. [19]. Failure is 
expected when, for a given mixed mode ratio G n / Gy, the calculated total energy release rate, Gy, exceeds the 
interlaminar fracture toughness, G c . Although several specimens have also been suggested for the measurement of 
the mode III interlaminar fracture toughness property [28-31], an interaction criterion incorporating the scissoring 
shear, however, has not yet been established. 
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Due to the lack of an existing fracture criterion for three-dimensional problems the in-plane energy release 
rates G//and Gm were combined to form the in-plane shearing component Gy=G//+G///, as suggested in [32]. Failure 
is now expected when, for a given mixed mode ratio G s /G r , the calculated total energy release rate, G/, exceeds the 
interlaminar Mode I/Mode II fracture toughness, G f . Therefore, for every computed mixed-mode ratio Gs /Gr along 
the delamination front as shown in Figure 8 the critical energy release rate G r , as shown in Figure 10, can be 
determined using the mixed-mode delamination criterion shown in Figure 9. Next, G? from Figure 7 was 
normalized by G r and was plotted as a function of specimen width as shown in Figure 1 1. Along the entire front the 
total energy release rate Gr calculated for damage onset load is much larger than the critical value G c (Gt/G c >\) 
which indicates unstable delamination growth. The critical load (P c ) which is required to advance a delamination 
from a pre-existing flaw is reached when Gr equals G r (Gj/G c -l). This type of analysis to obtain the critical load 
(Pc) has been performed earlier using two-dimensional models to asses the effect of assumed damage and location on 
the delamination onset predictions for skin-stiffener debonding [33]. 

The focus of the current study however, was the qualitative distribution of Gi/G c across the width of the 
specimen. The computed values indicate that in spite of peak total energy release rates at the free edge (z= 0.0 mm 
andz=25.4 mm) as shown in Figures 5 and 7, the delamination is not going to advance at the edges first, because 
the failure criterion is a strong function of the mode mixity. The distribution suggests that the delamination will 
progress asymmetrically in the interior first and that a slanted front will develop from an initially straight front 
which was assumed in the finite element model. An increased advance is to be expected with increasing z-values. 
This prediction is confirmed by X-ray photographs of a specimen taken during the tension fatigue testing as shown 
in Figure 12. The dark areas where the contrast fluid had penetrated into the damaged section of the specimen reveal 
the delaminations. In Figure 12(a) the gray intensity appears to increase towards the right (increasing z-values), 
which suggest an increased growth in this area. This observation is confirmed by Figure 12(b) where the 
delamination appears the have grown further on the right side of the specimen (increasing z-values). The correlation 
between prediction based on analysis and experiment appears encouraging, additional detailed modeling of the three- 
dimensional damage however appears necessary in order to fully understand the propagation mechanism. 

4. Application of a ShelI/3D Modeling Technique for the Analysis of Skin/Stiffener Specimens. The 

objective of this study was to demonstrate the usefulness of a shell/3D modeling technique for the investigation of 
delamination onset from an initial crack in the skin/stringer specimen. The shell/3D modeling technique developed 
earlier used a local three-dimensional solid finite element model only in the immediate vicinity of the delamination 
front, as shown in the example of a modeled Double Cantilever Beam (DCB) specimen in Figure 13(a) [17]. The 
approach combined the accuracy of the full three-dimensional solution with the computational efficiency of a plate 
or shell finite element model. Double Cantilever Beam (DCB), End Notched Flexure (ENF), and Single Leg 
Bending (SLB) specimens were analyzed in an earlier study using three-dimensional finite element models to obtain 
reference solutions [34]. Mixed mode strain energy release rate distributions were computed across the width of the 
specimens using the virtual crack closure technique. The analyses were repeated using the shell/3D technique to 
study the feasibility for pure mode I (DCB), mode II (ENF) and mixed mode I/II (SLB) cases [17]. For a local 
three-dimensional model, extending to a minimum of about three specimen thicknesses on either side of the 
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delamination front, the results were in good agreement with mixed mode strain energy release rates obtained from 
computations where the entire specimen had been modeled with solid elements as shown in Figure 13(b). The 
shell/3D modeling technique is believed to offer a great potential for reducing the model size, since only a relatively 
small section in the vicinity of the delamination front needs to be modeled with solid elements. The actual 
efficiency is described later. 

4.1. Shel!/3D Finite Element Model. The deformed three-dimensional model of the specimen with load and 
boundary conditions is shown in Figure 14. This full three-dimensional model was used to generate the reference 
solutions to which the results from the shell/3D simulations were compared [14]. The specimen was modeled with 
ABAQUS® solid twenty-noded hexahedral elements C3D20R with quadratic shape functions and a reduced 
integration scheme. To create the model, a two-dimensional mesh was extruded into ten elements across the width 
of the specimen as shown in Figure 14(a). A detail of the modeled delaminated region is shown in Figure 14(b). 
For short delamination lengths (a < 1.0 mm) crack opening mode I was apparent. For longer delaminations mode I 
ceased and the surfaces started to overdose as shown in the detail of Figure 14(c). This phenomenon had not been 
observed in the earlier results from two-dimensional analyses [33]. 

The deformed shell/3D model of the specimen is shown in Figure 15(a). The global section was modeled 
with ABAQUS® reduced integrated eight-noded quadrilateral shell elements S8R. The local three-dimensional 
model extended to about three specimen thicknesses on either side of the delamination front as shown in 
Figure 15(b). This configuration yielded good results during the initial feasibility study shown in Figure 13 [17]. 
The local three-dimensional section was modeled with solid C3D20R elements and was almost identical to the 
equivalent segment of the full three-dimensional model shown in Figure 14. A refined three-dimensional mesh was 
used in the critical area of the 45° skin ply where cracking was observed during the tests as shown in Figure 15(c). 
Outside the refined area, the mesh was modified to prevent the three-dimensional model from becoming excessively 
large. The skin plies were grouped into four layered elements with 45/-45/0, -45/45/90, 90/-45/45 and 0/45/-45 
respectively, thus taking advantage of the composite solid element option in ABAQUS® [21]. The fabric layers and 
the resin layer were grouped into five layered elements. In the transition regions several plies were modeled by one 
element with material properties smeared using the rule of mixtures [22]. This procedure did not ensure the full A- 
B-D contribution of the plies; however, it appeared suitable for the small transition regions to enforce a reasonable 
model size. The transition from the global shell element model to the local three-dimensional model in the vicinity 
of the delamination front was accomplished by using multi-point constraint options given by ABAQUS® to enforce 
appropriate translations and rotations at the shell-solid interface [21]. Using the shell/3D approach the total number 
of freedoms could be reduced by more than 12% compared to the full three-dimensional model of the skin/stringer 
specimen. 

The shell/3D model shown in Figure 15, however, was not fine enough in the vicinity of the free edges 
(2=0.0 mm and z=25.4 mm) to accurately capture the influence of the free edges on the distribution of the energy 
release rates across the width. Therefore, a shelI/3D model as shown in Figure 16 was created with twenty elements 
across the width of the specimen corresponding to the full three-dimensional model shown in Figure 4. The new 
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shell/3 D model was improved by a refined zone (0.74 mm, five elements) near the free edges (z=0.0 mm and 
z=25.4 mm) as shown in Figures 16(b) and (c). The efficiency is discussed in the summary. 

The Virtual C rack Closure Technique (VCCT) was used to calculate the mode I, II and III components of the 
strain energy release rates for the modeled delamination [18, 23]. For the entire investigation, the ABAQUS* 
geometric nonlinear analysis procedure was used. This was done in accordance with previous studies which were 
used as reference solutions [14]. 

4.2. Global Response. First, the global response of the specimens was computed at the mean quasi-static 
damage onset load determined from experiments. The load-displacement and the load-strain behavior computed from 
different FE models were compared to the corresponding experimental results. This global response was used to 
examine whether the FE models, the boundary conditions, the loads and the material properties used in the model 
yielded reasonable results. Strains were averaged from computed nodal point values over a length corresponding to 
the dimensions of the strain gages shown in Figure 2(a) [10]. 

The load versus displacement plot in Figure 1 7 shows that shell/3D model exhibits a stiffer behavior 
compared to results from the full three-dimensional model. The models also show a slightly stiffer response 
compared to the experiments (DCDT). This discrepancy may be explained by the fact that the material data used in 
the FE simulation originated from the literature. For a consistent simulation, material data should be taken from the 
batch of material that was used to manufacture the specimens. For comparison, results from 2D analysis were 
included in the plot of Figure 17 [14]. Results from the plane-strain analysis indicated a stiffer behavior as expected. 
This is caused by the constraints inherent to the plane-strain model, particularly in the ±45° plies. The plane-stress 
model which imposes the out of plane stresses to be zero and allows the displacement to be the free parameter 
exhibits - as expected - a more compliant behavior. The influence of the assumptions made in developing two- 
dimensional finite element models on skin-stiffener debonding specimens was studied in detail in [14] where 
geometrically nonlinear finite element analyses using two-dimensional plane-stress and plane-strain elements as well 
as three different generalized plane strain type approaches were performed. 

A comparison of measured strains at the surface of the flange (see Figure 2(a)) and computed results is shown 
in Figure 18. In Figure 19, measured strains at the surface of the top 45° skin ply near the flange tip (see 
Figure 2(a)) and computed surface strains were compared. For both locations, the strains calculated from the 
shelI/3D model are in good agreement with the results from the full three-dimensional finite element model. The 
models show a slightly stiffer response compared to the strains measured during the experiments. As before, results 
from 2D analysis were included in the plots of Figure 18 and 19 for comparison [14]. The plane-strain model 
showed a stiffer behavior yielding an upper bound while the plane-stress models were more compliant yielding a 
lower bound. 

4.3. Fractures Mechanics Analysis. The objective of this study was to demonstrate the usefulness of the 
shell/3D modeling technique for the investigation of delamination onset from an initial crack. A fracture mechanics 
approach was used and during a series of nonlinear finite element analyses, strain energy release rates were computed 
at each front location for the loads applied in the experiments. Results obtained form the shell/3D model are 
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presented in Figure 20 as a three-dimensional plot of the distribution of the total energy release rate, G 7 , across the 
width of the specimen. Values at the free edge (z=0.0 mm and z=25.4 mm) have been excluded from the plots as 
the model was not fine enough to accurately capture the influence of the free edges on the distribution of the energy 
release rates. Unlike the tension loading case it was observed that G r increases continuously along the length ( x - 
coordinate), except for the zones near the free edges (z=0.0 mm and z=25.4 mm) where a sharp decrease is observed 
caused by the overclosure of the interfaces after about 1 mm of delamination growth. The continuous increase of G T 
along the length indicates unstable delamination growth. Across the width (z-coordinate) the computed total energy 
release rate gradually increases with z before it drops off near the free edge (z=25.4 mm). 

A three-dimensional plot of the distribution of mixed mode ratio G s /Gr with delamination length is shown 
in Figure 21. Here Gy denotes the sum of the in-plane shearing components Gn+ Gnu and Gr denotes the total energy 
release rate G/+G//+G///, where G 7 is the opening mode. For two-dimensional analyses, where Gw= 0, this definition 
is equal to the commonly used definition of the mixed mode ratio, G n /G T . For three-dimensional analysis, which 
also yields results for the scissoring mode Gnu the modified definition of G s is introduced. As before, values at the 
free edge (z-0.0 mm and z=25.4 mm) have been excluded from the plots as the model was not fine enough to 
accurately capture the influence of the free edges on the distribution of the energy release rates. Along the length (x- 
coordinate) it was observed that the delamination initially starts with high shearing components, followed by a drop 
which is equivalent to an increase in opening mode I. For longer delaminations the shearing component remains 
constant. Unlike the tension loading case, the mixed mode ratio rises sharply at the edges for delaminations longer 
than 1 mm. This is caused by the overclosure of the delaminated surface near the edges of the model where the crack 
opening mode I disappears. A detailed study of this phenomenon would require contact analysis and preferably a 
refined model near the free edges. 

The total energy release rates along the centerline of the specimen (z-12.7 mm) obtained from three- 
dimensional and shell/3D analysis are plotted in Figure 22 for comparison. The values from the shell/3D analysis 
are within 6% of the results obtained from the full three-dimensional simulations. Improved results may be obtained 
by enlarging the section modeled with solid elements. The values from previous plain strain and plane stress 
analysis were included in Figure 22 for comparison [14]. Qualitatively all results follow the same trend. As before 
the values from plane strain and plane stress analysis form upper and lower bounds, except for very short 
delamination lengths. 

The variation of mixed mode ratio G s /G r with delamination length is shown in Figure 23. For short 
delamination length the mixed mode ratio yields high shearing components, followed by a drop which is equivalent 
to an increase in opening mode I. For longer delaminations the shearing component remains constant. The values 
from the shell/3D analysis are in good agreement with the results obtained from the full three-dimensional 
simulations. The values from previous plain strain and plane stress analysis were included in Figure 23 for 
comparison and follow the same trend [14]. The results from three-dimensional analysis in the center of the 
specimen width (z= 12.7 mm) show a higher shearing component compared to the results from two-dimensional 
analysis. 
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The distribution of the computed total energy release rate G T across the width (z-coordinate) of the specimen 
is plotted in Figure 24 for the longest delamination modeled (jc= 3 1 .2 mm). Across the width the computed total 
energy release rate gradually increases with z before it drops off near the free edge (z=25.4 mm). Results from a 
previous analysis obtained form the full three-dimensional model shown in Figure 14 have been included in the plot 
for comparison [14]. The G T values from the shell/3D analysis are 6% lower than the results obtained from the full 
three-dimensional simulations. The distribution of the mixed-mode ratio Gs IG T across the width (z-coordinate) of 
the specimen is plotted in Figure 25 for the same location (x=31.2 mm). Across the width the computed mixed 
mode ratio shows an almost constant value with increasing shearing components near the edges (z=0.0 mm and 
z=25.4 mm). 

For the longest delamination modeled (*=31.2 mm), the distribution of the computed total energy release 
rate G? across the width (z-coordinate) of the specimen is plotted in Figure 26. Results were obtained from the 
shelI/3D model with a refined mesh near the edges as shown in Figure 16. Across the width, the computed total 
energy release rate gradually increases before it first drops off near the free edges and then sharply increases in the 
zone of the mesh refinement in the immediate vicinity of the free edges (z- 0.0 mm and z=25.4 mm). As expected, 
the results obtained from the model with the refined mesh near the edges yield additional information about the 
steep increase in Gt which could not be captured with the less refined models (Figure 24). The distribution of the 
mixed-mode ratio Gs /G T across the width (z-coordinate) of the specimen is plotted in Figure 27 for the same 
location (jc= 3 1.2 mm). Across the width the computed mixed mode ratio is nearly constant and dominated by 
opening mode I but shows high shearing components near the edges which is caused by the overclosure of the 
delaminated surface near the edges where the crack opening mode I disappears. The current model provided a 
smoother distribution across the width and yielded additional information about the steep increase in mixed-mode 
ratio G t which could not be captured with the less refined model (Figure 25). 

A mixed-mode failure investigation as for the tension load case (see section 3.3) was not performed for the 
bending load case due to the lack of information regarding the size and shape of the delamination propagation. 
Further, a detailed study would also require additional contact analysis to avoid overclosure at the edges of the 
specimen during the simulations. 

4.4. Summary. Computed total strain energy release rates and mixed mode ratios obtained from shell/3D 
simulations were in good agreement with results obtained from full solid models. The concurrence of the results 
demonstrated the effectiveness of the shell/3D modeling technique for the investigation of delamination onset from 
an initial crack. However, for the current configuration a desired reduction in computation time could not be 
achieved. This is partially caused by the fine three-dimensional mesh in the vicinity of the delaminated region, 
which was not replaced by shell elements. Additionally, the transition from the global shell element model to the 
local three-dimensional model was accomplished by multi-point constraints to enforce appropriate translations and 
rotations at the shell-solid interface, which increased the computational effort. Nevertheless, for large built-up 
composite structures modeled with plate elements, the shell/3D modeling technique offers a great potential for 
reducing the model size, compared to a full three-dimensional model. In future applications the shell/3D modeling 
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technique will therefore progress towards structural components such as stringer stiffened panels as outlined in the 
overview chart of Figure 1. 

5. Concluding Remarks. The debonding of a skin/stringer specimen subjected to tension was studied using 
three-dimensional volume element modeling and computational fracture mechanics. A three-dimensional finite 
element model with a refined zone near the free edges was generated and delaminations of various lengths were 
discretely modeled at the locations where delaminations were observed during the experiments. The goal was to 
complement the results of an earlier investigation where the three-dimensional mesh in the vicinity of the free edges 
had not been fine enough to accurately capture the influence of the free edges on the distribution of the energy 
release rates across the width. 

Mixed mode strain energy release rates were calculated from finite element results using the virtual crack 
closure technique. The simulations revealed an increase in total energy release rate in the immediate vicinity of the 
free edges of the specimen which had not been captured by earlier models. The results also yielded an increase in 
shearing mode towards the free edge. Correlation of the computed mixed-mode strain energy release rates along the 
delamination front contour with a two-dimensional mixed-mode interlaminar fracture criterion became possible by 
introducing the in-plane shearing component Gs=G/^-Gui, and defining the mixed mode ratio as Gs t G T compared 
instead of the traditional two-dimensional definition of Gn t Gt . The approach suggested that in spite of peak total 
energy release rates at the free edge the delamination would not advance at the edges first because the failure criterion 
is a strong function of mode mixity. The qualitative prediction of the shape of the delamination front was confirmed 
by X-ray photographs of a specimen taken during testing. The good correlation between prediction based on analysis 
and experiment is encouraging future investigations. 

The application of the shell/3D modeling technique for the simulation of skin/stringer debond specimen 
subjected to three-point bending was also demonstrated. The study extended the application of this technique 
beyond the simulation of simple specimens - such as DCB, ENF and SLB specimens - where the delamination is 
located between unidirectional plies in the mid-plane of the specimen. The global structure was modeled with shell 
elements. A three-dimensional model, extending to about three specimen thicknesses on either side of the 
delamination front was used locally. Delaminations of various lengths were discretely modeled at the locations 
where delaminations were observed during previous experiments. 

Computed strains at the center of the flange were compared with results obtained from full three-dimensional 
finite element anaysis and experimental data. Results from both analyses and experiment were in good agreement. 
Calculated total strain energy release rates and mixed mode ratios obtained form shell/3D simulations also were in 
good agreement with results obtained from full solid models. The concurrence of the results demonstrated the 
effectiveness of the shell/3D modeling technique for the investigation of delamination onset from an initial crack. 
For comparison results from plane stress and plane strain models provided flange strains, as well as energy release 
rates, which formed an upper and lower bound of the results obtained from full three-dimensional or shell/3D 
simulations. 


ll 



The current study provided an additional verification step for the technique prior to its application to large, 
full-scale stringer stiffened panels because failure observed in the skin/stringer specimen is identical to the failure 
observed in the full-scale panels where the delamination is located between plies of different orientation. For large 
built-up composite structures modeled with plate elements, the shell/3D modeling technique offers a great potential 
for reducing the model size, compared to a full three-dimensional model, since only a relatively small section in the 
vicinity of the delamination front needs to be modeled with solid elements. 

Acknowledgements. The authors gratefully acknowledge Drs. T. Kevin O'Brien and Isabelle Paris for 
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TABLE 1 . 


Material Properties 


IM7/8552 Unidirectional Graphite/Epoxy Prepreg 

E\ i = 161.0 GPa 

£22 = 1 1.38 GPa 

£ 33 = 11.38 GPa 

v 12 = 0.32 

v 13 = 0.32 

v 23 = 0.45 

G\2 = 5.17GPa 

G ]3 = 5.17 GPa 

G 23 = 3.92 GPa 

IM7/8552 Graphite/Epoxy Plain Weave Fabric 

E\ i =71.7 GPa 

£22 = 71-7 GPa 

£ 33 = 10.3 GPa 

v 12 = 0.04 

v 13 = 0.35 

v 23 - 0.35 

G\2 = 4.48 GPa 

Gi3=4.14GPa 

G 23 = 4.14 GPa 

Grade 5 FM300 Adhesive 

£=1.72 GPa 

v = 0.3 

(assumed isotropic) 
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FIGURE 1 . Application of the shell/ 3 D modeling technique 


Q) 

■R ■o 

E | © 

o c o 
* 0=0 
CD CD CL 
<D -C X 

<0 Q) 


-r O 

o «•— 

CO "D 
CD p 


2 CO 
O) CL 


0) 

"O 

c 


CD o' 

2 co 

CO *- 
0 ) 3 

»- o 
u o 
c co 

Q) "O 

■e a) 

co co 

li 

I = 


CD 

| 

3 __ 

CL 0) 

Si 

£3 

o = 

3 3 

■o *♦- 

<D O 
■*—> 

C 

S £ 

O CO 
rE cl 
c £ 

D) § 
CO o 


CO 

c 

0) co 

C (D 

i o 

CD C 

c c 


c c 
.9 .9 

'*—> *—> L C 

U O O o 

3 3 3 -^ 

O O O o 

CD (D CD 3 

*o 

.o vO o 

o O ' o iZ 

CO CO CM O 



■D 


3 

E 

00 

o 

4 — 

■o 

c 

a> 

CD 

CD 

L_ 

h_ 

3 


O 

o 

c 

CO 


<D 

CD 

CD 

E 

O) 


a> 

T 3 

3 

CL 

E 

o 

c 

c 

o 



C 

CD 

C 

o 

Q. 

E 

o 

O 

"5 


o 

3 


CO 


1 


X 

CD 

0-1 

E 

o 

o 

*D 

<D 

CO 

03 

CD 

b 

c 


CO 

a> 

c 

o 

Gl 

3 

O 

o 


A 


E 

o 

■O 

CD 

a> 

i— 

o ^ 

<2 S ^ 


c 

c 8 

| E 

0 ) c 
c c 

c c 

o o 

V-* 

o o 

3 3 
“O -D 

a> <d 


i- o 

Tt 


CO 
<D CD 

£ E 

CD *+- 
CD -j 

■- o 

o -- 

•^3 C 

! 8 
*“ T 3 
0) 

o t_ 

N O 

Tt C 


c 

CD 

E 

0 
CD 

01 
CO 

c 

o 

cc 

N 

*k» 

a> 

o 

CO 

u. 

CO 

_c 

o 


15 


177.8 mm 


Tape skin 


1 .78 mm 


t 

A 



7 V 

Corner 3 " 42 mm Comer "I Strain Gages 


50.8 mm H 


Fabric flanae 



NOTE: Figure not to scale 


i 

3.86 mm 

A 


(a) Specimen configuration 



( b ) Tension test set-up 



(c) Three -point bending case 
FIGURE 2. Specimen configuration and test set-up . 
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{a) Deformed three-dimensional model of tension specimen with load and boundary' conditions 



( b) Detail of model of undamaged specimen 



(c) Detail of de formed model of damaged specimen 
FIGURE 4. Full three-dimensional model of skin/flange specimen 
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FIGURE 6. Computed mixed mode ratio G5/G7 for tension test. 
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FIGURE 9. Mixed-mode delamination criterion for IM7/8552. 
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FIGURE 10. Critical total energy release rate across the width of the specimen at a = 31.2 mm for tension test 
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FIGURE 12. Damage progression in skiti/st ringer debotid specimen 


23 



global shell element model 


of unfractured section 






global shell element model 


of upper arm 





















aMT 5 



detail of local 
3D FE model 
around 

delamination front 


■%r 


global shell element 
model of lower arm 


& jgl^ 


SI 


delamination front 
shell-solid transition interface 


(a) FE model and detail around de lamination front 
(c= 5 mm, d- 30 mm, e- 2 mm, /= 21.0 mm, 2h= 3.0 mm, B= 25.0 mm) 
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(b) Strain energy release rate distribution across the width 
FIGURE 13. Shell/3 D Finite element model of a DCB Specimen with [0]24 lavup. 
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(c) Detail of deformed model of damaged specimen 
FIGURE 1 5. Shell/3 D model of skin/flange specimen subjected to three-point bending 
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(b) Detail of undeformed local 3D model and shell to 3D interface region 



(c) Detail of deformed model of damaged specimen 


FIGURE 16. Shell/3 D model of skin flange specimen with fine mesh near the edge 
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FIGURE 17. Load-displacement plots for three-point bending tests 
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FIGURE 1 8. Flange strain-load plots for three-point bending tests 
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FIGURE 21. Computed mixed mode ratio G^/Gjjor three-point bending test. 
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